Counting statistics of heat transport in harmonic junctions — transient and steady 

states 
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We study the statistics of heat transferred in a given time interval tM, through a finite harmonic 
chain, caUed the center (C), which is connected with two heat baths, the left (L) and the right 
(-R), that are maintained at two different temperatures. The center atoms are driven by an external 
time-dependent force. We calculate the cumulant generating function (CGF) for the heat transferred 
out of the left lead, Ql, based on two-time measurement concept and using nonequilibrium Green's 
function (NEGF) method. The CGF can be concisely expressed in terms of Green's functions of 
the center and an argument-shifted self-energy of the lead. The expression of CGF is valid in both 
transient and steady state regimes. We consider three different initial conditions for the density 
operator and show numerically, for one-dimensional (ID) linear chains, how transient behavior 
differs from each other but finally approaches the same steady state, independent of the initial 
distributions. We also derive the CGF for the joint probability distribution P{Ql, Qr), and discuss 
the correlations between Ql and Qr. We calculate the total entropy flow to the reservoirs. In 
the steady state we explicitly show that the CGF obeys steady state fluctuation theorem (SSFT). 
Classical results are obtained by taking h ^ 0. The method is also applied to the counting of the 
electron number and electron energy, for which the associated self-energy is obtained from the usual 
lead self-energy by multiplying a phase or shifting the contour time, respectively. 

PACS numbers; 05.40.-a, 05.60.Gg, 05.70.Ln, 44.10.-l-i 



I. INTRODUCTION 



Nonequilibrium systems are common in nature because 
they are, in general, subject to thermal gradients, chem- 
ical potential gradients or may be triggered by time- 
dependent forces. Heat transport is one such example of 
nonequilibrium systems where the heat carriers could be 
electrons, phonons, magnons, etc. To study heat trans- 
port in phononic systems, one considers a finite junction 
part, which can be an insulator, connected with two heat 
baths that are maintained at different temperatures. In 
the past decade, the main focus was on the calculation of 
the steady state heat current or heat flux flowing through 
the junction part from the leads (l|-[To|. For diffusive sys- 
tems, the answer is given by Fourier's law [11 - 13] which 
is true only in the linear response regime, i.e., when the 
temperature difference between the baths is small. How- 
ever for harmonic or ballistic systenis, the heat current 
is given by a Landauer-like formula [1, 0, [l^l which was 
flrst derived for electronic transport. Landauer formula 
on the contrary to the Fourier's law is true for arbitrary 
temperature differences between the leads. No such ex- 
plicit expression for current is known for transient states. 
In recent times, several works [3, [l^ followed to answer 
what happens to current in the transient regime. This 
is an important question both from the theoretical and 
experimental points of view. 

Much attention has been given to phonon transport, 
in particular on thermal devices and on controlling heat 



flow [16[. With the advent of technology it is now pos- 
sible to study transport problems and observe a single 
mode of vibration in small systems with few degrees 
of freedom jl3|. These systems shows strong thermal 
fluctuations which play an important role because ther- 
mal fluctuations can lead to instantaneous heat transfer 
from colder to hotter lead. It is therefore necessary to 
talk about the statistical distribution of heat flux for 
these systems. In the electronic literature the distri- 
bution P{Ql) of the charge Ql, flowing from the left 
lead to the junction part, was answered by calculating 
the corresponding CGF, Z{£_) — (e*^'^^), and is given 
by the celebrated Levitov-Lesovik formula [T8l - [20j . This 
methodology is also known as the full counting statistics 
[2l| - l30| in the held of electronic transport. Experimen- 
tally the electron counting statistics has been measured 
in quantum-dot systems 31. , However few experi- 

ments have been done for phonons [s^ . In the phononic 
case Saito and Dhar [13] gave an explicit expression of 
the CGF. Ren et al. gave a result for two-level systems 
(35} . Full counting statistics of energy fluctuations in a 
driven quantum resonator is studied by Clerk [36| . The 
main focus in these papers was on the lon g-tim e limit 
and SSFT Usingthe NEGF method lH SI and 

two-time measurement |4l|, |4^, concept, Wang et 
al. [3^ gave an explicit expression for the CGF which is 
valid for both transient and steady state regimes. 

In this paper, we extend our previous work in Ref. [sol 
and derive the CGF in a more general scenario, i.e., in 
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the presence of both the temperature difference and the 
time-dependent driving force. We analyze the cumulants 
of heat Ql for three different initial conditions of the 
density operator and study the effects on both transient 
and steady state regimes. We also derive the CGF for 
the joint probability distribution of left and right lead 
heat P{Ql,Qr) which help us to obtain the correlations 
between Ql and Qr. By calculating CGF for P{Ql,Qr) 
we can immediately obtain the CGF for the total entropy 
that flows to the leads. We present analytical expressions 
of the CGF's in the steady state and discuss the SSFT. 
Our method can be easily generalized for multiple heat 
baths. 

The plan of the paper is as follows. We start in Sec. II 
by introducing our model. Then in Sec. Ill we define 
current and corresponding quantum heat operator fol- 
lowed by the definition of CGF for Ql using the two- 
time measurement concept, in Sec. IV. In Sec. V we 
derive the CGF using Fcynman's path integral method 
and in Sec. VI we use Feynman's diagrammatic tech- 
nique to derive the CGF. We discuss the steady state 
result and fluctuation theorems in Sec. VII. In Sec. VIII 
and IX we discuss how to calculate the CGF's numer- 
ically in transient regime and give numerical results for 
one-dimensional (ID) linear chain model, connected with 
Rubin heat baths, for three different initial conditions 
of the density operator. Then in Sec. X we obtain the 
CGF for joint probability distribution of heat transferred 
P{Ql,Qr) and discuss correlations and total entropy 
flow. In. Sec. XI we give the long-time limit expression 
for the driven part of the full CGF. In Sec. XII we discuss 
another definition of generating function due to Nazarov 
and discuss the corresponding long-time limit. We found 
that using this definition, the generating function does 
not obey the Gallavotti- Cohen (GC) fluctuation symme- 
try. We conclude with a short discussion in Sec. XIII. 
Few appendices give some details of technique nature. In 
particular, an electron system of a tight-binding model 
is treated using our method. 

II. THE MODEL 

Our model consists of a flnite harmonic junction part, 
which we denote by C, coupled to two heat baths, the 
left (L) and the right (R) , kept at two different tempera- 
tures Tl and Tr, respectively. To model the heat baths, 
we consider an infinite collection of coupled harmonic 
oscillators. We take the three systems to be decoupled 
initially and to be described by the Hamiltonians, 

na = lplPa + lulK''uo,, a = L,C,R, (1) 

for the left, right, and the finite central region. The leads 
are assumed to be semi-infinite. Masses are absorbed by 
defining u = \fmx. Ua and p„ are column vectors of 
coordinates and momenta. K" is the spring constant 
matrix of region a. Couplings of the center region with 



the leads are turned on cither adiabatically from time t = 
—CO, or switched on abruptly at t = 0. The interaction 
Hamiltonian takes the form 

Hint = ulV'^^uc + ulV''^uc. (2) 

For i > 0, an external time-dependent force is applied 
only to the center atoms, which is of the form 

Vc{t) = -f{t)uc, (3) 

where f{t) is the time-dependent force vector. The driv- 
ing force couples only with the position operators of the 
center. The force can be in the form of electromagnetic 
field. Coupling of this form helps us to obtain analytical 
solution for the CGF of heat flux. So the full Hamiltonian 
for t > (in the Schrodinger picture) is 

n{t) = H{o-)+Vc{t) = nc + nL+nR+Hint+Vcit). 

(4) 

In the next section we will define current operator and the 
corresponding heat operator based on this Hamiltonian. 

III. DEFINITION OF CURRENT AND HEAT 
OPERATORS 

It is possible to define the current operator X depend- 
ing on where we want to measure the current. Here we 

consider the current fiowing from the left lead to the cen- 
ter system and 1^ is defined (in Heisenberg picture) as 

(5) 

where T-LH{t) is the (time-dependent) Hamiltonian in the 
Heisenberg picture at time t. The corresponding heat 
operator can be written down as 

QlH) = flL{t') dt' = Hl{Q) - nl it), (6) 
Jo 

where Hl [ = ■Hz,(0)] is the Schrodinger operator of the 
free left lead and 

U^{t)=U{Q,t)ULU{t,0), (7) 

and lA{t,t') is the evolution operator corresponding to 
the full Hamiltonian T-L{t) and satisfies the Schrodinger 
equation 

^n^J^Ml = u{t)u{t,t'). (8) 

The formal solution of this equation is (assuming t > t') 

U{t,t')=Texpi^-^J^ H (t) dtj , (9) 

where T is the time-order operator where time increases 
from right to left. Also UHt,t') = U{t',t). Q of non- 
calligraphic font will be a classical variable. 

In the following section we derive the CGF based on 
this definition of heat operator and using two-time mea- 
surement scheme. 
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IV. DEFINITION OF THE GENERATING 
FUNCTION FOR HEAT OPERATOR 

Our primary interest here is to calculate the moments 
or cumulants of the heat energy transferred in a given 
time interval t^- Hence, it is advantageous to calculate 
the generating function instead of calculating moments 
directly. Since Ql is a quantum operator, there are sub- 
tleties as to how exactly the generating function should 
be defined. Naively we may use (e*^®^). But this defi- 
nition fails the fundamental requirement of positive def- 
initeness of the probability distribution. 

Here we will give two different definitions that are used 
to calculate the generating function for such problem. 
The first definition comes from the idea of two-time mea- 
surements and based on this concept the CGF can be 
written down as 



(10) 



which we will discuss in great detail in this section. 
The second definition of the CGF is 



i/2\ 



(11) 



where T is the anti-time order operator. The time (or 
anti-time) order is meant to apply to the integrand when 
the exponential is expanded and Qj^ is expressed as in- 
tegral over II as in Eq. This definition is used by 
Nazarov et al. [U mostly for the electronic transport 
case. In the last section we will show how this generat- 
ing function can be derived starting from Z(^) under a 
particular approximation and will also give explicit ex- 
pression for Zi(^) in the long-time limit. 

In the following we will discuss the idea of two-time 
measurement and derive the corresponding CGF Z{(^). 



A. Two-time measurement 

The heat operator in Eq. © depends on the left-lead 
Hamiltonian Hl at time and t. The concept of two- 
time measurement implies the measurement of a certain 
operator (in this case Hl) at two different times. Here 
the measurement is in the sense of quantum measurement 
of von Neumann \^ . 

Let us first assume that the full system is in a pure 
state I 'So) at t = 0. We want to do measurement of 
the energy associated with the operator Hl- According 
to quantum mechanics, the result of a measurement can 
only be an eigenvalue of the (Schrodinger) operator Hl 
and the wave function collapses into an eigenstate of "Hl . 
Let 



After the measurement at time t = 0, the wave function is 
proportional to H^ | ^I^o) if the result of the measurement is 
the energy a and the probability of such event happen is 
(^o|H^|^o). Let's propagate this state to time t and do a 
second measurement of the lead energy, finding that the 
result is b. The wave function now becomes proportional 
to IlbU{t,0)Ila l^'o)- The joint probability of getting a 
at time and b at time t is the norm (inner product) of 
the above (unnormalized) state. 

If the initial state is in a mixed state, we add up the 
initial probability classically, i.e., if 

p(0)=^«;fc|vl/g)(vl/g|, wk>0, ^«;fc = l, (13) 



the joint probability distribution of two-time measure- 
ment output is 



P{b,a) 



k 



Wk (*^^ \TlaU{0,t) UbUit, 0) H„ l^-g) 



Tr[Hap(0) UaUiO,t) UbUit^O)]. 



(14) 



By definition, we see that P{b, a) is a proper probability 
in the sense that P{b, a) > and b ^(^' = 1- Then 
the generating function for Ql = a — 6 is defined as 



a.b 



)p(6,a) 



J2 e*«^""'^ Tr [H, p(0) UaU{0, t) UbU{t, 0)] 



a,b 



where we define HJ 



p'(0)=^H,p(0)H,. 



(15) 



(16) 



which we call as the projected density matrix. 

If the initial state at < = is a product state i.e., 
p(0) = p(— oo) = PL ® PC ® PR, where the left, center 
and right density matrices are in equilibrium distribu- 
tions corresponding to the respective temperatures: pa — 
e-/3„H„/Xr[e-''°«°] for a = L,C,i? and /3„ = l/(/cBr„), 
then the projection operators H^ do not play any role 



and 



= Tr 



p(-oo) ■ 



Here we will derive the CGF for three different initial 
conditions: 

• Product initial state, i.e., p(— oo), which corre- 
sponds to sudden switch-on of the coupling between 
the leads and the center. 



a\4>a), Hq = \4>a){4>a\, 



(12) 



where H^ is the projector into the state \<j)a) satisfying 
H^ = Ha, and Ha = 1. We assume the eigenvalues are 
discrete (this is always so if the lattice system is finite). 



steady state as the initial state, i.e., p(0), which we 
can obtain, starting with the decoupled Hamiltoni- 
ans at t — —oo, switch on the couplings between 
the center region and the leads, adiabatically upto 
time t = 0. 
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• projected density matrix p'(0) considering p{0) as 
the steady state, i.e., taking the effects of measure- 
ments into account. 

In the following sections we will analytically show that 
the CGF's corresponding to different initial conditions 
reach the same steady state in the long-time limit and 
hence is independent of initial distributions. However 
for short time transient behavior depends significantly 
on initial conditions and also the measurements do play 
an important role. 



V. CALCULATION FOR 2(0 FOR INITIAL 
STATES p{0) AND p'(0) 

In this section we will give detail derivation for -H(^), 
using Feynman path-integral formalism, for two different 
initial density operators p{0) and p'(0). 



where h(x{t,t') is the modified evolution operator of an 
effective Hamiltonian given by 

= e"^^H(^)e~"^^ (21) 

where a; is a real parameter which in this case is — A 
and -C/2 - A. Finally Ux{t, t') is given by it > t') 

xe"'«^H(ti)H(t2) • • • H(t„)e^"«^ 
= Tei^v^~\j^n.{t')dt'Y (22) 

It is important to note that substituting A = in Z(^, A) 
gives us the initial density matrix p(0). 

Now we will give an explicit expression of the modified 
Hamiltonian which helps us to calculate the CGF 
using path integral. 



A. Removing the projection Ila at t = 



The projection by Ha at i = Eq. p6)) to the density 
matrix creates a problem for formulation in path inte- 
grals. We can remove it following Ref. liO by putting 
it into part of an evolution of "Hl, just like the factor 
associated with the generating function variable ^, with 
a price we have to pay, introducing another integration 
variable A. The key observation is that we can represent 
the projector by the Dirac 5 function 



Ha oc 6{a - Ul) 
dX 
2^' 



-iXia-nL) 



(17) 



For this to make sense, we assume the spectrum of the 
energy of T-Ll is continuous, which is valid if we take the 
large size limit first. Identifying Ha as S{a — Hl) with 
a continuous variable a introduces an constant propor- 
tional to the Dirac (5(0) to p'(0), since Ha is now nor- 
malized as HaHf, = S{a — 6)Ha. However, this constant 
can be easily fixed by the condition Z(0) = 1. So using 
Ha = S{a — Hl) will not cause difficulty. 

Substituting the Fourier integral representation into p' 
we obtain 



p'(0)cx y"daHap(0)Ha 



dX 
2^ 



e*^«XO)e 



(18) 
(19) 



Using the symmetric form of Z , Eq. ([T5|) , we have 
2(0 « / ^Tr{p{0)U^/,_^{0,t)U_^/,_^{t,0)} 



/dX 



(20) 



B. The expression for Tix 

The modified Hamiltonian is the central quantity for 
calculating CGF. It is the Heisenberg evolution of the 
full Hamiltonian Hit) (in Schrodinger picture) with re- 
spect to Hl- Since Ji^ commutes with every term H 
where ^.{t) = % + u^V^'^'uc, except the coupling term 
u we can write 

2"^^H(t)e-"^^ 
^H{t) + {uL{hx)^ULfv''^uc, (23) 



where ul^Hx) 
Heisenberg evolution to time t 
obtained explicitly as 



UL{hx) = cos(\/ KLhx)uL 



is the free left lead 
hx. UL{fix) can be 



:sm{y/K^hx)pL. (24) 



The matrix \^Kl is well-defined as the matrix Kl is posi- 
tive definite, ul and pl are the initial conditions att = 0. 
The final expression for T-Lx{t) is 



n^it) = nit) + [ulC{x) +pIS{x)\uc, (25) 



where 



rLC 



C{x) = {cOsihxyjKL) -I)V^ 

S{x) = {1/ ^/Kl)sm{hx^/KL)V^^ . 



(26) 
(27) 

The effective Hamiltonian now has two additional term 
with respect to the full H(t). The term u^C{x)uc is like 
the harmonic coupling term which modifies the coupling 
matrix V^'-^ . 

In the following we calculate the two parameter gener- 
ating function Z(^, A) using Eq. (PO)) . 
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C. Expression for Z{(,, A) 

The expression for A) can be written down on the 
contour as (see Fig. [T]) 



A) = Tr piO)T,e 



(28) 



where Tc is the contour-ordered operator which orders 
operators according to their contour time argument, ear- 
her contour time places an operator to the right. The 
contour function x{t) is defined as whenever t < or 
t > tM, and when Q < t < Im, i-e., within the mea- 
surement time interval, for upper branch of the contour 
x'^{t) = —£,/2 — A, and for lower branch x^{t) = — A. 

For the moment, let us forget about the other lead and 
concentrate only on the left lead and center. The effect 
of other lead simply modifies the self-energy of the leads 
additively, according to Feynman and Vernon (48| . Using 
Feynman path integral technique we can write 



Z(e,A)= / V[uc]V[uL]pi-^)e 



(i/h) fj^dT{Cc + CL+CLc) 

(29) 

Note that in Eq. the contour C is from to ti/ and 
back, while that in Eq. (I^H]) is on the Keldysh contour 
K, that is, from — oo to Im and back to take into account 
of adiabatic switch on, replacing p(0) by p(— oo). Their 
relation is 



p{0) ^U{0,-oo)p{~oo)U{~oo,0). 
We can identify the Lagrangian's as 



(30) 



C = Cl + Cc + C 



LC, 



Cr 



^LG = -ulSuc - ul{V^^ + C)uc. 



{K^-S^S)uc, 



(31) 



For notational simplicity, we have dropped the argument 
r. The vector or matrices /, C, and S are parametrically 
dependent on the contour time t. They are zero except 
on the interval Q < t < Im- / is the same on the upper 
and lower branches, while C and S take different values 
depending on x[t). 

Now the lead part can be integrated out by perform- 
ing Gaussian integral |48|. Since the coupling between 
the lead and the center is linear, it is plausible that the 
result will be a quadratic form in the exponential, i.e., 
another Gaussian. To find exactly what it is, we convert 
the path integral back to the interaction picture (with re- 
spect to T-Ll) operator form and evaluate the expression 
by standard perturbative expansion. The only difference 
is that now the coupling with the center involving both 
ul and ul- The result for the influence functional is 



FIG. 1. The complex-time contour in the Keldysh formalism. 
The path of the contour begins at time to, goes to time Im, 
and then goes back to time t — to- t and r' are complex-time 
variables along the contour, to = —oo and corresponds to 
Keldysh contour K and C, respectively. 



given by [i^ 

Il[u^{t)] EE J V[uL]pL{-^)ei 



f dT(CL+CLc) 



Tr 



e 21 



-Tee 



J dr J dT'u^(r)Il(^T,r')uc{T') 



(32) 



In the influence functional, the contour function uc (t) 
is not a dynamical variable but a parametric function. 
Vi(t) is the interaction picture operator with respect to 
the Hamiltonian H l and is given by 

V/(r) = pISuc + uliV'^'' + C)uc + ]^ulS'^Suc 

= uI{t + hx{T))V^'^uc + ]^u^S'^Suc. (33) 

The important influence functional self-energy on con- 
tour is given by 

n(r, r') = (r, r') + Ei(r, r') + S'^S S{t, r'), (34) 
(t,t') -I- S]l(t,t') = 1/^V(t + hx{r),T' + hx{r'))V'^^ 

^i:L{T + hx{T),T' + nx{T')), (35) 

where we obtain a shifted self-energy El (t + hx{T), t' + 
hx{T')'j which is the usual self-energy of the lead in con- 
tour time with arguments shifted by hx{T) and Hx^t'). 
We define the self-energy S;^ as the difference between 
the shifted self-energy and the usual one T,l{t,t'). T,^ 
turns out to be a central quantity for this problem as we 
will show that, the CGF Z can be concisely expressed in 
terms of the center Green's function Go and E^. 

Substituting the explicit expression for the influence 
functionals of both the left and right leads to the path 
integral expression given in Eq. (j29p . we have 



Z(e,A) - J 2?Mpc(-oo)e(^/'>)/''^^^/LM/i?M 
= / I?[uc]pc(-oo)e^'^"«, (36) 



(37) 



where the effective action is given by 



5e 



dr 



1 



1 



-\j dr j dT'ul{T){Y.{Ty) + ^t{Ty))uc[T'), 
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where S = + Si,;, taking into account the effect of 
both the leads. The S^S term in Il[uc] cancels exactly 
with the one in Cc- We can perform an integration by 
part on the ii^ term, assuming that the surface term does 
not matter (since it is at t = — oo), we can write the 
expression in a standard quadratic form 



SeS = ^ I dr I dT'ul{T)D{T,T')uc{T') 



+ j f{T)uc{T)dT. 

D{t, t') is the differential operator and is given by 

^(^' = -^|^'^(^' - ^^^^^^ 
-I](r,r')-Sf(r,r') 
= i?o(r,r')-Sf(r,r'). 



(38) 



(39) 



The above equation defines the Dyson equation on 
Keldysh contour. The generating function is obtained 
by doing another Gaussian integration and is of the fol- 
lowing form 



(40) 



(The meaning of the determinant will be explained in 
Appendix C). We define the Green's function G and Gq 
by DG = 1, and DqGq = 1, or more precisely 



D{t, r")G(r", T')dT" = IS{t, t'), (41) 



and similarly for Go . G can be written in terms of Go in 
the following Dyson equation form 

G{t,t')^Go{t,t') (42) 

dndT2Go{T, Tl)Sf (n, T2)G(t2, r'). 



We view the differential operator (integral operator) 
D and as matrices that are indexed by space j and 
contour time t. / is a column vector. The exponential 
factor term can also be written as a trace, f^D~^f = 
TiQT){Gff^). We can fix the proportionality constant 
by noting that Z(^ = 0, A = 0) = 1. Since when ^ = 0, 
A = 0, we have a; = and thus (r, r') = T,l{t + x,t' + 
x') — Si(T, r') = 0, so D = Dq. The properly normalized 
CGF is 

Z(C, A) = det{D^^D)-^/^e-^^^^''f. (43) 

We don't need to do anything for the exponential factor 
because of the following reason. We note 



fGof^ J J dTdT'f{TfGo{T,T')f{T') (44) 

= Y.I I ^dtcj'dt'f{tfGf{t,t')f{t'). 



Since the driven force / does not depend on the branch 
indices, i.e., f~^(t) = f~{t), we can take the summation 
inside and obtain 



(Ta'G""' = G^ + G*, - G> - G< = 0. 



(45) 



Finally making use of the formulas for operators or 
matrices det(M) = e'^''^"*^, and ln(l — y) = ~J2T=i V 
we can write the CGF in terms of S;^ for the projected 
initial condition case as, 

lnZ(0 = lim lnZ(^, A) 



lim 

X—yoc 



oo 



^ 2n 

n—1 



(GoSf) 



(i,r)(G'oSf GoSf ) + 



2h 



/^GoS^Go/ 



(46) 



This expression for CGF is valid for any transient time tM 
present in the self-energy and is the starting point for 
the calculation in transient regime. The notation Tr^- 
means trace both in space index j and contour time r (see 
Appendix C). In order to obtain Z{£^) from Z(^, A) we 
have to take the limit A — > oo because Z{^, A) approaches 
a constant as |A| — >■ cxd and hence the value of the integral 
is dominated by the value at infinity. Since S^(t, t') = 
for ^ = we have the correct normalization Z{0) — 1. 

Similarly, for the steady state initial condition p{0) the 
CGF is given by 

lnZ(0 = linilnZ(e,A) (47) 

A-s-O 

The difference in this two cases is in the matrix E;^. 

Similar relations also exist if we want to calculate the 
CGF for right lead heat operator Qr. In this case one 
has to do two-time measurement on the right lead corre- 
sponding to the Hamiltonian Hr. The final formula for 
the CGF remains the same except S;^ should be replaced 

by ^i- 

Now in order to calculate the cumulants {{Qa)) with 
a = L, R we need to go to the real time using Langreth's 
rule |44| . In this case, it is more convenient to work 
with a Keldysh rotation (see Appendix C) for the contour 
ordered functions while keeping Tr(Ai3G ■ ■ ■ D) invariant. 
The effect of the Keldysh rotation is to change any given 
matrix T)'^'^ {t,t'), with a,(j' — zL for branch indices, to. 



V 



pi? pa 



(48) 



~ 2 + V* - V< ~ V> , V< - 2?* +V* -V> 

In this case we define the quantities I?"", P", , and V'^ 
as above. In particular, ^ 2?< + 2?>, as one usually 
might thought it is. 
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Using the above definition for the center Green's func- 
tion Go we get 



Go 



Gr r^K 
^0 

Gg 



(49) 



The component is due to the standard relation 
among Green's functions. But the K components are 
not zero for and Go, as we will compute later. 

It is useful computationally to work in Fourier space 
even if there is no time translational invariance. We de- 
fine the two-frequency Fourier transform by 



dt 



+ 00 



dt'A{t,t')e 



l\ i{uit+ui't') 



Since Go is time-translationally invariant then, 
Go [w, uj'] = 27r5 (cj + w')Go [uj] , 



(50) 



(51) 



is "diagonal" , where the single argument Fourier trans- 
form is similarly defined. 



A[lo] 



+ 00 



(52) 



(The expressions for different components of Go[a;] and 
S [lo] are given in Appendix A) . Using Gq [uj] , we can save 
one integration due to the d function, and have 



InZ(e) 



2 ^ 



Tr 

2h '^^''^ 



dn 1 - GoMSf 



(53) 



where Go[a;]S;^[w, w'] is viewed 

as a matrix indexed by lo and lo'. The trace is per- 
formed on the frequency as well as the usual space and 
branch components. (The meaning of trace in frequency 
domain is discussed in Appendix C). T is given by 



:F[lo,lo'] 



2/M/M^ 





(54) 



In the next section we derive the CGF for the product 
initial condition using Feynman diagrammatic technique. 
Because of the special form of the initial density matrix 
the calculation for the CGF simplifies greatly in this case. 



VI. PRODUCT STATE p(-oo) AS INITIAL 
STATE 



In this section, we derive the CGF starting with a prod- 
uct initial state, i.e., the density matrix at time t = is 
given by p(—oo) = Pc®Pl®Pr- Since this density matrix 
commutes with the projection operator IIq, the initial 
projection does not play any role in this case. Working 
in the interaction picture with respect to the decoupled 



Hamiltonian T-L{—oo) = X^i^j; the interaction part of 
the Hamiltonian on the contour G = [0, Im] is 

VJW = -f{T)uc{T)+UR{T)V'''''uc{T) 

+ ul{t + hx{T))V^^uc{T). (55) 

In the last term for u^, the argument is shifted by hx 
where x+{t) = -^/2, x~ {t) = i/2 for < i < Im- 

The density matrix remains unaffected by the transfor- 
mation to the interaction picture, because it commutes 
with H{—oo). The CGF can now be written as 



(56) 



Expanding the exponential, we generate various terms of 
product of Ua ■ These terms can be decomposed in pairs 
according to Wick's theorem [3]. Since the system is 
decoupled, each type of u comes in an even number of 
times for a non-vanishing contributions because (uc) — 
0, (ucul) — and we know 



{TcUa{T)Ua'{T') ) p{-r^) ^ Sa,a' 9a{T, t' ) . (57) 



We use Feynman diagrammatic technique to sum the 
series, since V/ contains only two-point couplings, the 
graphs are all ring type. The combinatorial factors can 
be worked out as l/(2rt) for a ring containing n vertices. 
We use a very general theorem which says In Z contains 
only connected graphs, and the disconnected graphs can- 
cel exactly when we take the logarithm. The final result 
can be expressed as 



lnZ(a = ~iTr,-,ln 



1 - 



(58) 



where 



= SL(r + x, t' + a;') + Sj?,(t, r') = S + , (59) 

and S is the total self-energy due to both the leads. 
G(r, t') obeys the following Dyson's equation 



G{T,T')=gc{T,T') 



(60) 



dTidT2gc{r, ri)I]'°*(ri, r2)G(T2, r'). 

The above expression for CGF can be written down 
more explicitly, by getting rid of the vacuum diagrams. 
Let us define a new type of Dyson's equation 

Go{r,r')^gc{r,r') (61) 

dri dr2 5c (t, Ti ) S (n , r2 ) Go (r2 , t' ) , 



where gc is the contour ordered Green's function of the 
isolated center. (The Green's functions for an isolated 
single harmonic oscillator is given is appendix A). This 
expression looks formally the same as before except that 
Go satisfies a Dyson equation defined on the contour from 
to tM and back, while G is defined on the Keldysh 
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contour from — oo to tM- Using this definition we can 
write 



l-5cS'°*-l-.gc(S 



= (l-gcS)(l-GoSf). (62) 

The two factors above are in matrix (and contour time) 
multipHcation. Using the relation between trace and de- 
terminant, Indet(M) = TrlnM, and the fact, det{AB) = 
det(A) det(i?), we find that the two terms give two fac- 
tors for Z, and the factor due to 1 ~ gc^ is exactly 1. 
We have then 



lnZ(0 = -^Tr,,.ln 



1 - GoSf 



(63) 



where the G{t,t') can now be expressed in terms of 
Go{t,t') as 



G- 



Gn 



yA 



(64) 



which is similar in form to Eq. (|42|) . 

The expression for In Z(^) is consistent with the earlier 
result, Eq. (1461) . in the long-time limit. So we can con- 
clude that the long-time limit is the same independent of 
the initial distributions. 

To compute the cumulants ((Q")), we need to take 
derivative with respect to ^ n-times to \nZ. Note that 
the shifted self-energy for < t < tjv/ is (for all three 
initial conditions) 



S^(t,t') 



yt 



S<(i,0 
E>(t,t') 



0, 
0, 



(65) 



We note Sf (^ = 0) = 0. The derivatives at f = can be 
obtained as 



9e 



4=0 



= (-;i)"s<^^"^(i- 



(i-O, 



(66) 



where the superscript (n) means derivatives with respect 
to the argument of the function n times. In the following 
sections we first show the explicit expression of the CGF 
in the long-time limit and then discuss the steady state 
fluctuation theorem. 



to all the Green's functions and self-energy can be per- 
formed (where the translational invariance is restored). 
Applying the convolution theorem to the trace formula 
in Eq. (j63p . we find that there is one more time integral 
left with integrand independent of t. This last one can 
be set from — tA//2 to tA//2, obtaining an overall factor 
of Im and we have 



Tr(,-,)(AB---i?) = iM 



2^ 



-Tr 



A{u:)B{uj) ■ ■ ■ D{uj) 



(67) 

In the long-time limit, the shift given to the argument 
in depends on the branches, and the two arguments 
(t, t') becomes t — t' and we have 



^A 



{t,t') ^ Y^l" {t+x^-t'-x" ) 
S$(<) = I]>(i + fte)-S^(i). 



L 



{t~t'), (68) 
(69) 



Fourier transforming the greater and lesser self-energy, 
we get 



I]>H(e 



1) 



(70) 
(71) 



We note that E;^ is supposed to depend on both ^ and A. 
However in the long-time limit, the A dependence drops 
out which makes the steady state result independent of 
the initial distribution. 

Finally, we can express the generating function as 



InZ(C) 



duj 

— Trln 
47r 



1 



47r 



-Tr 



G[uj]P[uj, 



(72) 



where G[a;] is obtained by solving the Dyson equation 
in frequency domain and in the long-time obeys time- 
translational invariance. So the full CGF can be written 
as the sum of contributions due to driving force and due 
to temperature difference between the leads, i.e., 



lnZ(0 =lnZ^(e)+lnZ''(e). 



(73) 



In the following and subsequent sections we discuss about 
Z^{^) and we will return to Z'^{S_) in Sec. XI. 

In order to obtain the explicit expression for \\iZ''{^) 
we need to compute the matrix product 



VII. LONG-TIME LIMIT AND STEADY STATE 
FLUCTUATION THEOREM 

For the long-time limit calculation we can use either 
Eq. ([53]) or Eq. (|63| . For convenience of taking the large 
time limit, i.e., large, we prefer to set interval to 
(—tM I'^T't'M 1'^)- In this way, when Im oo, the in- 
terval becomes the full domain and Fourier transforms 



GoHsfH 



K 



Gr y^iS 

Gg 



'{a + h) b — a 



(74) 



To simplify the expression, we rewrite the term Trln(l — 
M) as a determinant and use the formula (assuming A 
to be an invertible matrix) 



det 



A B 
C D 



detiA)det{D -GA-^B) (75) 
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to reduce the dimensions of the determinant matrix by 
half. The steady state sohition for Z^{(^) is given by 

Indetjl-GSriGSr^ 



47r 



2)/l/7?]}. (76) 

with fa = l/(e'3°'*'^ - 1), a = L,R, the Bose-Einstein 
distribution function and Ta[oj] = — S^[a;]). If 

we consider the full system as a one-dimensional linear 
chain, then because of the special form of Fq matrices 
(only one entry of the F matrices are non-zero) it can be 
easily shown that 

det[/ - (GsriGgF^)s(e)] - 1 - rHs(e) (77) 

where S(^) is any arbitrary function of ^ and T[w] — 
Tr(GorLGgFfl) is known as the transmission function 
and is given by the Caroh formula The generating 

function 2*(^) in the steady state obeys the following 
symmetry 



(78) 



where A = /3r — /3l is known as thermodynamic affin- 
ity. This relation is also known as Gallavotti-Cohen (GC) 
symmetry [37] . The immediate consequence of this sym- 
metry is that the probability distribution for heat trans- 
ferred Ql which is given by the Fourier transform of the 
CGF, i.e., P{Ql) = ^ JZ, d^2iO e-'^'?- obeys the fol- 
lowing relation in the large tM limit. 



(79) 



This relation is known as the steady state fluctuation 
theorem and was first derived by Saito and Dhar [34^ in 
the phononic case. This theorem quantifies the ratio of 
positive and negative heat flux and second law violation. 

The cumulants ((Q")) can be obtained by taking 
derivative of In (^) with respect to and setting ^ = 0. 
The first cumulant is given by 



(Ml 

tM 



/CO 1 



(80) 



which is known as the Landauer-like formula in ther- 
mal transport. Similarly the second cumulant ((Q^)) = 
(Q^) ~ {Q)^: which describes the fluctuation of the heat 
transferred, can be written as [3^ [sollsij. 

{huj) |r [uj) {fL - /r) 



tM 



4tt 



+ T{UJ) {fL + fR + 2fLfR)}. 



(81) 



Our formalism can be easily generalized for multiple heat 
baths and for N leads connected with the center G, we 
can generalize the above formula as 



lnZ^(0 - -tM 



J ^ lndct{/-^GJiriGf;F„[(e^«'"--l) 



In the following section we will discuss the how to numeri- 
cally calculate the CGF in the transient case for projected 
density matrix p'(0). We also discuss about solving the 
Dyson equation given in Eq. (j6ip . 



VIII. TRANSIENT REGION 

The central quantity to calculate the CGF numerically 
is the shifted self-energy T,f^ which is given by 

j:i{T,T') ^^L{T + hx{T),T' + hx{T')) ~^l{t,t'). (83) 

Here r is a contour variable which runs over Keldysh con- 
tour K — (—00,00) and back, for the initial conditions 
p(0) and p'(0), whereas for p(— oo), r runs over the con- 
tour G = [Oj^m] (see Fig. 1). The contour function x{t) 
is whenever t < or t > tM, and ior < t < tM, 
x+(i) = — ^/2 — A, and (t) ^ £_/2 — X. Depending on 
the values of i, t', and A (A — >■ and A — >■ oo corresponds 
to steady state initial state and projected initial state, 
respectively) will have different functional form. If 
< t,t' < tM then Ef's are given by Eq. (gS]). This 
is the region which dominates in the long-time limit and 
gives steady state result. If both t and t' lies outside the 
measurement time, i.e., t,t' <0 or t,t' > tM then E;^ is 
zero. 



The main computational task for a numerical evalu- 
ation of the cumulants is to compute the matrix series 

- ln(l - M) ^ M + ^AP -\ . It can be seen due to 

the nature of T,^ that for the product initial state, exact 
n terms upto M" is required for the n-th culumants, as 
the infinite series terminates due to S;^(^ = 0) = 0. Nu- 
merically, we also observed for the projected state p'(0), 
exactly 3n terms is required (although we don't have a 
proof) if calculation is performed in time domain. 

The computation can be performed in time as well as in 
the frequency domain. However for projected and steady 
state initial condition since Go[i^] is time translational 
invariant it is advantageous to work in the frequency do- 
main. But for the product state there is no such prefer- 
ence as Go in Eq. (j6ip is not time translational invariant 
and one has to solve it numerically. 

In the following we first discuss how to calculate 
T,f^[ijj,uj'] for projected initial state, defined in Eq. ([55)) 
and then we will discuss how to solve the Dyson equation 
for the product initial condition case, given in Eq. (j6ip . 



fL+ (e-'«'"''^-l)/™ + (e^«''"+e-''«''"-2)/L/„] } 



A. calculation of E^(a;, Lj') 

To calculate T,^{uj,uj') for projected initial state p'{0) 

we define two types of theta functions 9i{t,t') and 
02{t,t'). 9i{t,t') is non-zero when 

^ 0<t<iM, and t' < or t'>tM, (84) 
or 

(82) 0<t'<tM, and t<0 or t>tM, (85) 
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and 92{t,t') is non-zero only in the regime where < 
t,t' < Im ■ For the regions where 6i{t,t') is non-zero 
the expression for T,^ after taking the limit A oo is, 
(assuming all correlation functions decays to zero as i — > 
±00) 



{t,t' 



it-t'). 



(86) 



So using theta functions we may write T.^{t,t') in the 
full t,t' domain as 

i]^(i, t') = ~ei{t, t')s< {t - 1') + 02(i, t') X 
[^<{t-t'-hO-^<{t-t')\ 

Y.> {t, t') = -0i {t, t')Y.> {t - t') + e2{t, t') X 

[Y>{t-t' + hi)-^>{t^t')\ (87) 
By doing Fourier transform it can be easily shown that 



= - / — (71 



3r 



W-Wc, w'-l-Wc] S^*(lUc) 



and 



/OO 7 

/OO 1 
-co 

where r/ — ±1. The positive sign is for and negative 
sign for S^. 

The theta functions are now given by 

6'i(a;Q, Wb) = j{uja)-g{^h) + f{ujb)-g{i^a), 

92{u:a,LOb) = .f{L0a):f{0Jb), (90) 

where 

/H = 



(91) 

ILO + € lUJ — e 

with e — > 0"*" . The theta functions are of immense impor- 
tance which carries all information about the measure- 
ment time tM- 

In the limit Im — >• oo, the region Q < t,t' < Im dom- 
inates and corresponding theta function, i.e., 02(w,w') 
reduces to 

92{UJ — UJctUj' + UJc) ~ 5{u) — LJc)5{uj' + Wc), (92) 

and is responsible for obtaining the steady state result. 

To calculate all the cumulants we only need to take 
derivative of Yia{uj^lu') with respect to since Gq does 
not have any ^ dependence. Also E'^ has ^ dependence 
only for Q <t,t' <tM and hence the derivatives are given 
by 



(?7?IClIc)"6'2 [ui — OJc,io' + OJc 



E>'<(tJc)e""=''^. (93) 
Here n refers to the order of the derivative. 



B. Dyson equation on contour C 

Let us now discuss about solving the Dyson's equation 
for Go given in Eq. (I6ip for product initial state p{—oo). 
In order to compute the matrix Go{t,t') we have to cal- 
culate two components Gq and Gq which are written in 
the integral form by applying Langreth's rule [3, [EO] 



GS(t,i')=5&(i-0 



(94) 



+ / dh dt2gh{t-ti)^''{ti-t2)Gl{t2,t'), 
Jo Jo 



and 



G^{t,t')=g^{t-t') 

+ dh dt2ghit-ti)E'^ih~t2)G^it2,t') 
Jo Jo 

+ dh dt2ghit-ti)J:''ih-t2)G^{t2,t') 
Jo Jo 



(95) 



+ / dh dt2gBit-h)^''ih~t2)G^„it2,t'). 
Jo Jo 

Note that the argument for center Green's function gc 
and lead self-energy E are written as time difference t — t' 
because they are Green's functions for isolated center 
part and leads respectively and hence are calculated at 
equilibrium. The analytical expressions for E and gc are 
known in frequency domain and are given in Appendix A. 
To determine their time-dependence we numerically cal- 
culate their inverse Fourier transforms using trapezoidal 
rule [52|. Then in order to solve above equations for any 



Im we discretize the time variable into N total intervals 
of incremental length At = Im/N and thus converting 
the integral into a sum. After discretization, the above 
equations can be written in the matrix form which are 
indexed by space j and discrete time t, as 



Gil ryr\^K /^a 
— ^0^ ^0 



(I + Gl^^Tgcil + ^^G^o)- (96) 



So Gq can be obtained by doing an inverse of the matrix 
(/ — ^pE*") and then multiplying by g^. Gq in this case 
also obeys time-translational invariance, so it can also be 
obtained by direct inverse Fourier transform. Gg can be 
obtained by taking transpose of Gq. Once Gq and Gq 
are obtained we use the second equation to calculate Gq^ 
which is simply multiplying matrices. 

Similarly E^ in Eq. (j65p are obtained by doing inverse 
Fourier transforms of the lead self-energy. We follow the 
same steps independently to calculate the cumulants for 



IX. NUMERICAL RESULTS 

We now present some numerical results. In Fig. 2 and 
3, we show the results for first four cumulants for both 
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FIG. 2. (Color online) The cumulants {{Ql)) and {{Qr)) for 
n=l, 2, 3, and 4 for one-dimensional linear chain connected 
with Rubin baths, for the projected initial state p'(0). The 
black and red curves corresponds to ^^nd {{Q'r)) re- 

spectively. The temperatures of the left and the right lead 
are 310 K and 290 K, respectively. The center (C) consists of 
one particle. 




20 40 60 



FIG. 3. (Color online) Same as in Fig. 2 except for product 
initial state p(— oo). The temperatures of the left, the center 
and the right lead are 310 K, 300 K and 290 K, respectively. 




Im(10 s) 



FIG. 4. (Color online) Same as in Fig. 2 except for steady 
state initial state p(0). 



Ql and Qn (measurement is on the right lead) for ID 
hnear chain connected with Rubin baths, starting with 
the projected initial state p' (0) and product state p(— oo) 
respectively. Rubin baths |53l Is^ j mean in our case a 
uniform linear chain with spring constant K and a small 
onsitc Kq for all the atoms. Only one atom is considered 
as the center. The atoms of the left and right side of the 
center are considered baths. We use K = 1 eV/(uA2) 
and the onsite potential Kq — 0.1 eV/(uA^) in all our 
calculations. First of all, cumulants greater than two are 
nonzero, which confirms that the distribution for P{Ql) 
or P{Qr) is not Gaussian. The generic features are al- 
most the same in both the cases. However the fluctu- 
ations are larger for the product initial state p(— oo) as 
this state corresponds to the sudden switch on of the 
couplings between the leads and the center and hence 
the state is far away from the correct steady state dis- 
tribution. On the contrary, for the initial state p'(0) the 
fluctuations are relatively small. For p'(0) due to the 
effect of the measurement, at starting time energy goes 
into the leads, which is quite surprising. But for p(— oo) 
although initial measurement do not play any role, en- 
ergy still goes into the leads. This can also be shown 
analytically (see Appendix B). At the starting time the 
behavior of both Ql and Qn are very similar and can be 
understood since both the left and right leads are identi- 
cal and the effect of temperature difference is not present. 
However at longer times the odd cumulants starts differ- 
ing and finally grows linearly with time tjv/ and agrees 
with the corresponding long-time predictions. 

In Fig. 4 we show the results for the steady state initial 
condition, i.e., p(0) which can be obtained by mapping 
the projection operators as identity operator, i.e., taking 
the limit A — > 0. So in this case measurement effect is 
ignored and the dynamics starts with the actual steady 
state for the full system. The first cumulant increases Hn- 
early from the starting, {Q) = tl and the slope gives the 
correct prediction with the Landauer-like formula. How- 
ever, high order cumulants still have transient behavior. 
In this case the whole system achieve steady state much 
faster compared with the other two cases. 



X. CORRELATION BETWEEN LEFT AND 
RIGHT LEAD HEAT 

A. Product initial state 

In this section, we derive the CGF for the joint proba- 
bility distribution P{Ql, Qr) for the product initial state 
p(— oo). In order to calculate the CGF we need to mea- 
sure both T-Ll and T-Ln at time and at time Im- Since the 
Hamiltonians for the left and the right lead commute at 
the same instance of time i.e., \Hl,'Hi{\ — 0, such type 
of measurements are allowed in quantum mechanics and 
also Nelson's theorem [s^ gurentee's that P{Ql,Qr) is 
a well-defined probability distribution. The immediate 
consequence of deriving such CGF is that, the correla- 
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tions between the left and the right lead heat can be 
obtained and it is also possible to calculate the CGF for 
total entropy flow (defined below) to the reservoirs. To 
calculate the CGF we need two counting fields and £,r 
and the CGF in this case can be written down as ^] 

(97) 



where the average is defined as 



• )' = ^n^nfp(o)n^nf. 



and are the projectors onto the eigenstates of 
and Hji with eigenvalues a and c respectively, cor- 
responding to the measurements at t = 0. Here we will 
consider only the product state p(— oo), then initial pro- 
jections and Trf- do not play any role. We can proceed 
as before and finally the CGF can be written down as 



E ^Tro,.) 
fe=i 



i.e., in this case we need to shift the contour-time ar- 
guments for both left and right lead self-energies. In the 
long-time limit Z{£^l,^ji) becomes a function of difference 
of counting field and i.e., — fi?- The explicit 
expression for the CGF in the long-time limit is 

In Z(eL - ^r) = -tM J ^ Indctj/ - GSriGgPfl 

[(e«(«^-««)''"_l)/^ + (e-^(«^-««"'"'-l)/fl 

+ {e'(ir.-iR)nu.^e-^ia-iR)nu_2)f^fj^]y (loO) 

where Go obeys the same type of Dyson equation as in 
Eq. (|6ip . This CGF in the steady state obeys the same 
type of GC fluctuation symmetry, which in this case is 
given by 



(101) 



Now performing Fourier transform of the CGF, the 
joint probability distribution is given by P{Ql,Qr) = 
P{Ql) 5{Ql + Qr)- The appearance of the delta func- 
tion is a consequence of the energy conservation in the 
steady state, i.e., Il = —Ir- In the steady state knowing 
probability distribution either for Ql or Qji is sufficient 
to know the joint probability distribution. 

The cumulants can be obtained from the CGF by 
taking derivatives with respect to both and ^fl, 
i.e.,((Q£Q^)) = 9"+"MnZ/a(i^L)"9(i^fl)'", substitut- 
ing — — 0. In the steady state the cumulants obey 

{{qiqr)) = (-i)'"((Qr")) = The 

first cumulant give us the left and right lead correlation 
{{QlQr}) = {QlQr)-{Ql){Qr) and in the steady state 
is equal to —{{Q\))- 

In Fig. 5 we plot the first three cumulants for one di- 
mensional linear chain connected with Rubin bath where 
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FIG. 5. (Color online) First three cumulants of the cor- 
relations between left and right lead heat flux for one di- 
mensional linear chain connected with Rubin baths, starting 
with product initial state p{—oo). The left graph corresponds 
to {{QlQr)) and the right graph corresponds to cumulants 
{{QIQr)) (Black curve) and {{QrQl)) (Red curve). The left, 
center and right lead temperatures are 310 K, 290 K and 300 
K respectively. The center (C) consists of one particle. 




(10 s) 



FIG. 6. The cumulants of entropy production ((cr")) for n=l, 
2, 3, 4 for one dimension linear chain connected with Rubin 
baths, for product initial state p(— oo). The left, center and 
right lead temperatures are 510 K, 400 K, and 290 K respec- 
tively. The center (C) consists of one particle. 



the center consists of only one atom. Initially the cu- 
mulant {{QlQr)) is positively correlated as both Ql 
and Qli are negative, however in the longer time since 
Ql — —Qr the correlation becomes negative. We also 
give plots for {{Q\Qr)) (black online) and {{Q^^Ql)) 
(red online) which are in the long-time limit negative 
and positively correlated respectively and match with the 
long-time predictions. 



B. Entropy flow to the reservoir 

From the two parameter {^l,£,r) CGF one can also 
obtain the total entropy that flows into the leads. The 
total entropy flow to the reservoirs can be defined as (sgI . 



(t = -PlQl- PrQr- 



(102) 
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In order to calculate this CGF we just make the substi- 
tutions £_L — >■ —f^LfJ- and —Pr^J- in Eq. (|^ . In the 
long-time limit the expression for entropy-production is 
similar to ln2^(^L,^ij) with — replaced by A and 
becomes an explicit function of thermodynamic affinity 
Pr — Pl [l^- The CGF in this case satisfies the following 
symmetry 



(103) 



In Fig. 6 we give results for the first four cumulants of 
the entropy fiow. All cumulants are positive and in the 
long-time limit give correct predictions. 



XI. LONG-TIME RESULT FOR In 

In this section we derive the explicit expression for the 
long-time limit of the CGF lnZ''(^) which is given by 
(Eq.ESl) 



duj 
An 



Tr[GMJ"[w,-c 



(104) 



where G[uj] obeys the Dyson equation given in Eq. (|42l) . 
It is possible to write down G[uj] in terms of Go and 

as G[uj] = (/ — Go^f) ^Go[ll>]. This equation can be 
solved analytically. Next we assume that the product of 
/(<) and f{t') is a time-translationally invariant function, 
i.e., f{t)f'^{t') = F{t ~ t') in order to get rid of t + 
t' dependence term. In the Fourier domain this means 
= 2'kF[uj\5{uj + uj'). So from Eq. ^ the 
matrix element F12 is given by ^[w,— cj]i2 cx 5{Q)F[ijj\. 
We write (5(0) = tul'^'^- Using these results the CGF 
can be expressed as 



(105) 

where a and h are defined in Eq. ([70)) and Eq. ([7T|) . Using 
the expressions for the self-energy the CGF reduces to 



In 2-^(0 = 



Tr 



GSHr4^]GgMFH 



with 



(106) 



^(0 = (e-'^'^'^-l) + /L(e'«^'" -2), (107) 



and 



AA(e) = det[/ - (GSriGgrfl){(/ie»«'"--l) -f /fl 

{e-'^^'^-l) + (e''«'""^+e-*«'*"-2)/L/i?j.(108) 



It is important to note that /C(^) depends only on left lead 
temperature and satisfies the symmetry /C(^) = /C(— ^ — 
iPh)- So we can immediately write Z'''{—i(3L) = 1 and 
this relation is completely independent of the information 



about the right lead. If we consider the two leads at the 
same temperature (/3l — /3fl — /3), this form of symmetry 
is then closely related to the Jarzynski equality (JE) [5^, 
[soj and Z'^{—iP) = 1 is one special form of JE. However 
since M{C) does not satisfy this particular symmetry of ^ 
at thermal equilibrium (it obeys the GC summery when 
the leads are at different temperatures) and the CGF 
\aZ'^{S_) doesn't satisfy any such symmetry relation and 
hence JE is not satisfied. This does not violate JE as 
our definition of Z'^{^) is different from the one used to 
derive JE. 

Let us now come back to the general scenario with 
leads at different temperatures and give the explicit ex- 
pression of first and second cumulant by taking derivative 
of lnZ''(^) with respect to i^. 

The first cumulant or moment is given by (60j 



(109) 



where we define S[uj] as the transmission function for the 
driven case and is given by 



5H =Tr[GSriGgi^]. 



(110) 



From the expression of S[uj] it is clear that the average 
energy current due to driven force is independent of h and 
since it contains Gq" and F^, which are independent of 
temperature we can conclude that the energy current is 
independent of the temperature of the heat baths in the 
ballistic transport case. However the second cumulant 
and similarly the higher ones do depend on temperature 
of the baths. The second cumulant can be written as 



{{QD) 

tM 



2r(w)(/L-/fl) 



(111) 



Similarly all the higher cumulants can be obtained from 
the CGF and we can conclude that the distribution 
P{Qd) is not Gaussian. 



A. Classical limit 

In this section we will give the classical limit of the 
steady state expression for the CGF In Z^{^) and lnZ''(^) 
given in Eq. ^ and Eq. (11(15)) . 

First of all we note that retarded and advanced Green's 
functions, i.e., Gq and Gg are similar both for quantum 
and classical case, so they stay the same when fi, — 0. 
We know that in the classical limit fa — > ^f^"" ^^'^ ^^^^ 

e^^ — 1 + ix + — K • • • , where x = ^fvjj. Using this we 
obtain from Eq. ([76]) the classical limit of Z^{S). 



\tiZ: 



duj Indet 



/ — {G^qTlGqTr) 



kBTLkBTRi^i^ + A) 



(112) 
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This result reproduces that of Ref. [42| . In the classical 
case also the CGF obeys the GC symmetry, i.e., it re- 
mains invariant under the transformation — > — — yl. 

Let us now get the classical limit for \nZ'^{£^) using 
Eq. (jl06l) . Following above relations the function /C(^) 
in the limit ?i — > reduces to 



(113) 



The transmission function S[uj] stays the same as it is 
independent of temperature and h. So in the classical 
limit lnZ'*(^) reduces to 



+ 



(114) 



where 



det 



/ - (GSFiGgrj^) kBTL ksTR 



(115) 



Here we can easily see that Z'^^—iPl) = 1. 

We can also derive the fluctuation dissipation theorem 
from Eq. (|llll) if we assume the leads are at the same 
temperature, i.e., Pl — Pb, — P then we can write the 
second cumulant {{Q^)) as 



HQD) 
t]\i 



(M'5M(1 + 2/l) 



(116) 



In the high-temperature limit using /l — > '^fy^^ and we 



obtain 



Pj 



-{Qd). 



(117) 



In the next section we discuss Nazarov's generating 
function and give long-time limit expression. 



XII. NAZAROV'S DEFINITION OF 
GENERATING FUNCTION 

In this section we will derive another definition of CGF 
given by Eq. ([TT|). starting from the CGF, derived using 
two-time measurement concept, i.e., Eq. (jlOp . Eq. (fTTj) 
can be obtained from Eq. (jlOp in the small ^ approxima- 
tion as follows. In the small ^ approximation the modified 
Hamiltonian given in Eq. (|25p takes the following form 



(118) 



because lim^_j.oC(a;) = and lim^_^o'5(x) — hxV^^ . Xl 
is defined in Eq. ([5]). So the modified unitary operator 
becomes 

U^{t, 0) = Te-i m(t)+nxir.mdi _ (119) 



We can consider TixIl (0) as the interaction Hamiltonian 
and write the full unitary operator lAx as a product of 
two unitary operators as following 



U,{t,Q)=U{t,Q)Ui{t,Q), 



(120) 



where 



Z^(t,0) =re-* 

with lL{t') = W{t\0)lL{0)U{t',0) is the current oper- 
ator in the Heisenberg picture. It is important to note 
that U is the usual unitary operator which evolves with 
the full Hamiltonian H{t) in Eq. ([9]) and has no ^ depen- 
dence. 

If we use product state p(— oo) as the initial state the 
CGF is given by 

Z(0 = Tr[p(-(»)Z^^/2(0,i)Z^_5/2(i,0)]. (122) 

In the small ^ approximation and using the expressions 
for Uj. we can write the CGF as 

Zi(0 = ImZ(e) = Tr[p(-oo)Z^//2(0,i)Z^!e/2(i,0)], 

(123) 

where we use the property of unitary operator, i.e., 
W {t,0)U{t,0) = 1. Finally using the definition of heat 
operator Ql given in Eq. ^ and the CGF can be written 
down as 



Zi(0 = ^Te'«eL(t)/2j.g^?CL(t)/2 



(124) 



which is the same as in Eq. (ITTT) . 

In the following we will give the long-time limit expres- 
sion for this CGF. 

In order to calculate the CGF, it is important to go to 
the interaction picture with respect to the Hamiltonian 
"Ho = 'Hl+'Hc+'Hr, as we know how to calculate Green's 
functions for operators which evolves with 'Ho and treat 
the rest part as the interaction Vx — 'Hint + fi.xXi(0). So 
the CGF on contour C — [0,^*/] can be written as 



(125) 



where V^(t) is now given by 

Vi(T) = ul{T)V'^^uc{r) + ul{r)V^^uc{r) 

+nx{T)pL{T)V''^uc{T), (126) 

where pl = u^. The time-dependence r is coming from 
the free evolution with respect to T-lo- x{t) has the similar 
meaning as before, i.e., on the upper branch of the con- 
tour x+(t) — — ^/2 and on the lower branch x^ {t) = £_/2. 
Now using the same idea as before, we expand the series, 
use Wick's theorem and finally the CGF can be expressed 
as 



lnZ(0 = -^Tr,-.ln 



1-GoS 



(127) 
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Here Go is the same as before and is given by Eq. (1611) . 
However the shifted self-energy S;^ in this case is different 
and is given by (in contour-time argument) 

^iir, t') = hx{T) T,p^uL (t, t') + hx{T') T^u^^pi^ (r, r') 

+ft2 2;(r)x(T')Sp,pJr,r'). (128) 

The notation E^b(t, r') means 

^ab{t,t') = (-^)F^^(r,A(T)i3^(T'))F^^. (129) 



The average here is with respect to equihbrium distri- 
bution of the left lead. It is possible to express the 
correlation functions such as Ep^^j^ (r, r') in terms of 
the S„^,„^(t,t') = El(t,t') correlations. Sp^„^(t,t') 
and S„^p^(r, r') is simply related with T,l{t,t') by the 
contour-time derivative whereas for Sp^^p^ (t, r') the ex- 
pression is 



drdr' 



+ 5(T,T')Ei. (130) 



Where E£ 



yCLyLC _ ^qw in the frequency domain 
A 



different components of E;^ takes the following form 



yt 
^A 






f^'eyi 

4 


yi 
^A 




?72c2, ,2 


r^'eyi 

4 


y< 

^A 




— [inEuj 

V ^ 4 




^A 




= — infijj 

V 4 





(131) 



Finally using the relations between the self-energy (see 
Appendix A), in the long-time limit the CGF can be 
written down as, 

InZi(C) = -tM / ^ In [1 - {i^n^nM Ul - Jr) 



(rH(l + 2/L)(l + 2/i^)-GgE 

(132) 



where J'(^^,^'*) is given by 



1 {i^huj)^ 

4 4 



l(m 
4 4 



GSEiGgEi. (133) 



This CGF does not obey the GC fluctuation symmetry. 
However it gives the correct first and second cumulant 
as it should because the definition of first and second 



cumulant turn out to be the same for both the generating 
functions Z{£^) and -Zi(^) and is given by 



- (Q) = 



d\nZ{C) 91nZi(0 



= [ dti ( dt2{lL[ti)lL{t2)) - 

Jo Jo 



dm m) JO 



dh{lL{h)), 



dhilLih)) .(134) 



Expressions for higher cumulants are different for the two 
generating functions and hence the final expressions for 
the CGF's are completely different from each other. 



XIII. CONCLUSION 

In summary, we present an elegant way of deriving 
the CGF for heat Ql,r transferred from the leads to 
the center for driven linear systems using the two-time 
measurement concept and with the help of the NEGF 
technique. The CGF is written in terms of the Green's 
function of the center and the self-energy E^ of the leads. 
The counting of the energy is related to the shifting in 
time for the self-energy. This expression is valid in both 
transient and steady state regimes, where the informa- 



tion about the measurement time tM is contained in E 
The form of the expression, -(l/2)Trln(l-GoEf), is the 
same whether we use a product initial state or a projected 
initial state, except that the meaning of the Green's func- 
tion has to be adjusted accordingly. We consider three 
different initial conditions and show numerically for ID 
linear chains connected with Rubin baths, that transient 
behaviors significantly differs from each other but even- 
tually leads to the same steady state distribution in the 
long-time limit. We give explicit expressions of the CGF 
in the steady state invoking the symmetry of transla- 
tional invariance in time. The CGF obeys the GC sym- 
metry. We also give the steady state expression for the 
CGF in the presence of time-dependent driving forces. 
We obtain a two parameter CGF which is useful for cal- 
culating the correlations between heat flux and also the 
total entropy which flows to the leads. Our calculations 
can be easily generalized to arbitrary dimensions with 
any number of heat baths. We will show in the appendix 
that our method can be extended for the electronic cal- 
culations where we derive the CGF for a tight-binding 
model. It will be interesting to derive the CGF by tak- 
ing magnetic field contribution into the Hamiltonian and 
also to study the cumulants in the presence of nonlin- 
ear interactions such as phonon-phonon interactions or 
electron phonon interactions. 
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APPENDIX 



Expressions for different type of Green's 
functions 



Here we give the exphcit expressions for the center 
Green's function Go[u>] in the steady state, for a har- 
monic system which is connected with the leads. These 
formulas are required to derive the analytical form of the 
CGF given in Eq. (fTS)) . For the basic definitions of dif- 
ferent types of Green's functions we refer to Ref. il3- 

The retarded Green's function G'Q[a;] is given by 



(135) 



Here rj is an infinitesimal positive number which is re- 
quired to satisfy the condition of causality i.e.,GQ{t) — 
for t < 0. The advanced Green's function is Gq[w] = 

[GSM]^ The Keldysh G reen's function Gq[uj] can be 
obtained by solving the corresponding Dyson equation, 
Eq. (j6ip . and is given by 



where Ef + E| and Ef = E< + E> with a = 

L,R. Alternatively, Gg = Gq+Gq. Another important 
identity is 

GSM - GSM = -zGSM(rLM + r/?M)GSM, (i37) 

where r„M = «(5^aH - ^aH), and a = L,R. The 
self-energy for the leads are given by 

E<M = ./aM(s:;M-ssH), 

E>M = (1 + /oM)(s:,M - KH- (138) 

where /q[w] — l/(e'^°''"° — l) is the Bose distribution 
function. 

Explicit expressions for Gq[u;] and E£[a;] can be ob- 
tained for ID homogeneous linear chain, with inter par- 
ticle force constant K and onsite spring constant Kq and 
which is divided into three parts: the center, the left and 
the right. The classical equation of motion for the atoms 
in all three regions is 



ilj = Kuj-i + (^—2K — Ka)uj + Kuj^i, 



(139) 



where the index j runs over all the atoms in the full 
system. 

The retarded Green's function Gq[w] can be obtained 
by solving Q [(w -|- irj)^ — K]Gq = /, where matrix K 
which is infinite in both directions and is 2K + Kq on the 
diagonals and —K on the first off-diagonals. The solution 
is translationally invariant in space index and is given by 



g; 



X\3-k\ 



OJk I 



(140) 



with A = -^±^^/W^^AK^ and Vl = {uj+iT]f-2K- 
Kq, choosing between plus and minus sign by |A| < 1. 

The surface Green's function (7£[a;] can be similarly 
obtained in frequency domain and is given in terms of 
the self-energy E£[w] — —KX. Since in equilibrium only 
one Green's function is independent, knowing E£[w] is 
sufhcient to obtain all other Green's functions. 

Here we also give the expressions for Green's functions 
gc in time and frequency domain for an isolated single 
harmonic oscillator with frequency ujq (we have omitted 
the subscript G in gc) [Ml, [gj] 



fit) = ~9it) 



sincjQ^ 

Wo 

1 



9 M = 



2 ' 



-ITT 



g<M [S{Lj + ojo){l + f)+S{oj-ojo)f] ,(141) 



where / = /(wq) 



Other components can 



be obtained by exploiting the symmetry between the 
Green's functions such as g°'{—t) = g'^{t) for t > hence 
^'"[a;] = g"[— oj]. The greater component is related with 
the lesser component via g^{t) = g'^{—t) which in the 



(136) frequency domain satisfy ^^[a;] = g^[- 



B. Current at short time for product initial state 

p(-oo) 

Using the definition of current operator given in Eq. ([5]) 
the energy current flowing from the left lead to the center 
is (here we assume that there is no driving force f{t)) 

(Mt)) = -(^^) = 'j:{[nut),n]), (142) 

where the average is with respect to /?(— oo). If t is small 
we can expand T-Lhit) in a Taylor series and is given by 
Hl(0 =HL(0)+^HL(0) + ••• 
Now since [/9(— oo),'Hl(0)] = 0, then it immediately 
follows that ([Hl(0),?^]) = by using the cyclic prop- 
erty of trace. So in linear order of t the current is given 

by 



{idt)) ^t'-{[nL{Q),n]) ^ ~tl{[plv^^uc,n]). 



(143) 

The only term of full T-L that will contribute to the is 

Now using the relation that [pl,ml] ~ —ih, for one- 
dimensional linear chain we can write 

(140) = -tK^iu^)') = -tK'^^[fciuo) + l). 

,(144) 

where uf is the first particle in the center which is con- 
nected with the first particle of the left lead with force 
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constant K. Now since the average is with respect to 
/9(— oo),((uf )^) can be easily computed. Here fci^o) is 
the Bose distribution function of the particle with char- 
acteristic frequency wo- So we can see that for short 
time the current is negative, i.e, it goes into the lead. 
It is now easy to see that similar expression should also 
hold for {Xji{t)). The negative sign in currents means 
that the energy flows into the leads initially irrespect to 
the temperature of the leads. This is consistent with the 
numerical results obtained by Cuansing et al. [H, [l^ . 



C. Convolution, trace, and determinant on 
Keldysh contour 

Here we discuss the meaning of convolution, trace and 
determinant on the Keldysh contour which we used to 
derive the CGF's for heat flux. We define the convolution 
on contour in the following way. 

AB---D^ j dT2--- j dTnAj,j^{Ti,T2) 

32 ,31. ,3n 

), (145) 

From the convolution we define trace by substituting 
Tn+i = Ti, jn+1 = ji and integrate also over ti, sum 
over ji i.e., 

Tr,- , {AB---D) = J dTi y dT2 • • • J dTn (146) 

Tt, [A(ti, T2)B(r2, ra) • • • i?(T„, n)] , 

Changing from contour to real-time integration from — cxi 
to -|-oo, i.e., using J dr = I '^^^ have 

Tij^riAB ■■■D)= ^ j dti j dt2--- j dtn (147) 

,^2 ^<yn 

Tr, [A^i'^^ {tut2)^2B"'"'' (^2, is) • • • a„i?""'^"+i (t„, ii)] . 

Let us absorb the extra a into the definition of branch 
components, i.e., define 

A„„, ^ aA'^''', or A = a,_A, (148) 

where A is viewed as 2 x 2 block matrix with the usual 
+, — component. 



(149) 



(150) 



Then we can do a rotation, where the rotation matrix is 
given by 



A = 

and Gz is defined as 



( A++ 


A+- ^ 




(A^ 


A<\ 


[a-+ 


A- , 




[a> 


A* j 



1 

-1 ; ' 



then it can be easily seen that 

TTj^r{AB---D)^ J dti j dt2--- j dtr,TTj[A{ti,t2) 
B{t2M) ■ ■ ■ D{tnM)\ 

= TTt^j,a{AB---D). (151) 



(152) 



and we define for any matrix A, the rotated matrix as 

A = O^a^AO = O'^AO. (153) 

This is known as Keldysh rotation. The effect of Keldysh 
rotation is given in Eq. Since this is an orthogonal 

transformation the trace remains invariant and hence we 
can write 

Trtj- ,(^B ■■■D)^ Trt,,- ,(iB ■■■!)). (154) 

If we now go to the frequency domain using the definition 
of two-time Fourier transform given in Eq. (j50p then we 
can compute the trace in frequency domain as 

A[uji, -a;2]-B[a;2, -t^s] • ■ ■ -D[w„, - wi]} 
= Tr,-,,^(iB---^). (155) 

The last line above define what we mean by trace over 
frequency domain given in Eq. (j53p . Unlike trace in time 
domain, the second argument of the each of the variables 
need a minus sign. 

Let us now define what do we mean by 1 on contour. 
In the sense of convolution we define 1 as 



AID ^ AD 



(156) 



which means 



j dnj dr2A{T,n)I5{Tl,T2)D{T2,r')^ J dTlA{T,Tl)D{n 

(157) 

Note that 6{t, t') in the real time has the following form 
5"^"' {t,t')^a5^^a-5{t-t'). (158) 
The inverse on the contour is defined as 

j dTiAiT,Ti)B{Ti,T') = IS{t,t'), (159) 

where the identity matrix I takes care about the space 
index. Similar to the above we go to the real time and 
multiply the above equation with the branch index a and 
we can write, 

J dtiA{t, ti)B{h,t') = 15{t - t'). (160) 

where 

5{t - t') ^ aS"-"' {t, t') = (7'^5„,^,5{t - t') 

= 5,^,,5{t~t') (161) 



18 



If we now discretize the time and write 5{ti, tii) = Si^i' /At 
with At = It, — tit I then we have 



AB = /. 



(162) 



with A = A At and similarly for other matrices. 

With similar notions we can now write different types 
of Dyson's equation given in Eq. (|42l6ip as following. In 
contour time we have 



Go(T,r')=5c(T,T') (163) 

+ dTldT2gc(T,Tl)i;(Tl,T2)G'o(T2,T'), 

In real time following the above arguments we write 
Go(t,t') =5c(i,t') (164) 
+ dt^dt-2 gc{t,ti)f:{tut2)Goit2,t'), 

After Keldysh rotation we can write 

Goit,t') ^ gc{t,t') (165) 

+ J J dtidt2 gc{t,ti)t{tut2)Go{t2,t'). 

Finally in the discretize time t we write 

Go ^ 9c + ~9c^Go, (166) 

which is a matrix equation. Similar equations can also 
be written down for Eq. I^^ . 

Now we define determinant via the relation det(j4) — 
exp(TrlnA), i.e, the determinant is defined in terms of 
trace. In order for In A to be defined we have to assume a 
Taylor expansion. For example we can define ln(l+M) = 

M - ]VP/2 + H where 1 means %'(5(t,t') in 

contour space. 



D. A quick derivation of the Levitov-Lesovik 
formula for electrons using NEGF 

The generating function for the non-interacting elec- 
trons was first derived by Levitov and Lesovik [3 
using Landauer type of wave scattering approach. Klich 
[23 | and Schonhammer ^2§\ re-derived the formula us- 
ing a trace and determinant relation to reduce the prob- 
lem from many-body problem to a single particle Hilbert 
space problem. Esposito et al. gave an approach using 
the superoperator nonequilibrium Green's function (40| . 
A more rigorous treatment is given in Ref. l63l . 

Our method for calculating CGF can be easily ex- 
tended for the electron case. Here we will derive the 
CGF for the joint probability distribution for particle 
and energy without time-dependent driving force. The 
Hamiltonian of the whole system can be written as (us- 
ing tight-binding model) 



a=L,C,R 



clh"c^+ J2 (clK"^cc+h.c.) (167) 



where Cq is a column vector consisting of all the anni- 
hilation operator of region a. c]j, is a row vector of the 
corresponding creating operators, is the single par- 
ticle Hamiltonian matrix. V"'^ has similar meaning as 
in the phonon Hamihonian and V^^ = (Vf^^Y ■ 
We are interested in calculating the generating func- 
tion corresponding to the particle operator A/^ and en- 
ergy operator Hl oi the left- lead where Hl = c\h^CL 



and A/l — c]^cl |6^. One can easily generalize the for- 



mula for right lead also as we did in the phonon case. For 
electrons JVl and Hl can be measured simultaneously 
because they commute, i.e., [Hl,J^l] = 0. In order to 
calculate the CGF we introduce two counting fields 
and for particle and energy respectively. Here we will 
consider the product initial state (with fixed tempera- 
tures and chemical potentials for the leads) and derive 
the long-time result. 

Similar to the phonon case we can write the CGF as 



2(?e,Cp) 



(168) 



where superscript H means the operators are in the 
Heisenberg picture at time t. In terms of modified Hamil- 
tonian the CGF can be expressed as 



where 



V 2 ' 2 



,)(0,t)Z^(_ 



.^)(i,0)). 



-im 



(169) 



(170) 



with X = and y = ^p/2 and U{t,0) - 
is the modified Hamiltonian which evolves with both H l 
and Al and is given by 



"Ha; . 



He 



— ixl-LL—'iy-^L 



^Ul+Hc +'Hr + {e'ycl{hx)V^^cc + h.c.) 
+ (4K^^cc + h.c.), (171) 
where we have used the fact that 

e""^CL(0)e-"«^ = CL{hx), 



JyJ^L 



(172) 

So the evolution with Hl and A/l is to shift the time- 
argument and produce a phase for cl,c^ respectively. 
Next we go to the interaction picture of the modified 



Hamiltonian 'Hx,y with respect to Hq 



a=L,R 



and the CGF then can be written on the contour running 
from to tM and back as, 

Z(ee,ep) - Tr[p(-oo)T,e-Hrf-vi.„(r)j^ (^73) 

where y{T) is written in contour time. 

Kyi^) = (e''c[(r + fe)K^^cc(r) + h.c.) + 

(cfl(r)tT//c^cc(r)+h.c.). (174) 
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Now we can expand the exponential in the generating 
function and use Feynman diagrams to sum the series 
and finally the CGF can be shown to be 



ln-2(^e,^p) = Trj.^ln 



1 



Ge sr^A 



(175) 



where we define the shifted self-energy for the electron 
case as 

,(T,r') = e^(^'(^')~^(^))Ei,e(T+fi:E,r'+fi,^')-SL,e(r,r'). 

(176) 

The counting of the electron number is associated with 
factor of a phase, while the counting of the energy is 
related to translation in time. Note that the CGF does 
not have the characteristic 1/2 pre- factor as compared 
to the phonon case because c and c' are independent 
variables. In the long-time limit following the same steps 
as we did for phonons, the CGF can be written down as 
(after doing Keldysh rotation) 

— Trln {l-Gl{E)tl,{E)). 

(177) 

In the energy E domain different components of the 



shifted self-energy are 



e 



\{E) = 0, 



^«=^)-i)i]>(£;). 



Finally the CGF can be simplified as 
\nZ = tM lndet{/ + Gj;rLGgrfl 



+ (e- 



l)/fl 



-2)/l^ 



(178) 

"-1)/l 
}. (179) 



where a = £,p+£,eE and /q, is the Fermi distribution. Note 
the difference of the signs in the CGF as compared to 
the phonons. If we replace ahy [E — fJ-L)^, the resulting 
formula is for the counting of the heat Ql = Hl — ^J■L■^L 
transferred, where /x^ is the chemical potential of the left 
lead. The CGF obeys the following fluctuation symmetry 

m 

z{^,,^p) ^ z{-^,+i{fiR-pL),-^p-i{-(3Rm-hf^L))- 

(180) 
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